home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / sqrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  9.1 KB  |  410 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    sqrt   double precision square root
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    sqrt
  29.  *    machine independent routines
  30.  *    math libraries
  31.  *
  32.  *  DESCRIPTION
  33.  *
  34.  *    Returns double precision square root of double precision
  35.  *    floating point argument.
  36.  *
  37.  *  USAGE
  38.  *
  39.  *    double sqrt (x)
  40.  *    double x;
  41.  *
  42.  *  REFERENCES
  43.  *
  44.  *    Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
  45.  *
  46.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  47.  *    1968, pp. 89-96.
  48.  *
  49.  *  RESTRICTIONS
  50.  *
  51.  *    The relative error is 10**(-30.1) after three applications
  52.  *    of Heron's iteration for the square root.
  53.  *
  54.  *    However, this assumes exact arithmetic in the iterations
  55.  *    and initial approximation.  Additional errors may occur
  56.  *    due to truncation, rounding, or machine precision limits.
  57.  *    
  58.  *  PROGRAMMER
  59.  *
  60.  *    Fred Fish
  61.  *
  62.  *  INTERNALS
  63.  *
  64.  *    Computes square root by:
  65.  *
  66.  *      (1)    Range reduction of argument to [0.5,1.0]
  67.  *        by application of identity:
  68.  *
  69.  *        sqrt(x)  =  2**(k/2) * sqrt(x * 2**(-k))
  70.  *
  71.  *        k is the exponent when x is written as
  72.  *        a mantissa times a power of 2 (m * 2**k).
  73.  *        It is assumed that the mantissa is
  74.  *        already normalized (0.5 =< m < 1.0).
  75.  *
  76.  *      (2)    An approximation to sqrt(m) is obtained
  77.  *        from:
  78.  *
  79.  *        u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
  80.  *
  81.  *        P0 = 0.594604482
  82.  *        P1 = 2.54164041
  83.  *        Q0 = 2.13725758
  84.  *        Q1 = 1.0
  85.  *
  86.  *        (coefficients from HART table #350 pg 193)
  87.  *
  88.  *      (3)    Three applications of Heron's iteration are
  89.  *        performed using:
  90.  *
  91.  *        y[n+1] = 0.5 * (y[n] + (m/y[n]))
  92.  *
  93.  *        where y[0] = u = approx sqrt(m)
  94.  *
  95.  *      (4)    If the value of k was odd then y is either
  96.  *        multiplied by the square root of two or
  97.  *        divided by the square root of two for k positive
  98.  *        or negative respectively.  This rescales y
  99.  *        by multiplying by 2**frac(k/2).
  100.  *
  101.  *      (5)    Finally, y is rescaled by int(k/2) which
  102.  *        is equivalent to multiplication by 2**int(k/2).
  103.  *
  104.  *        The result of steps 4 and 5 is that the value
  105.  *        of y between 0.5 and 1.0 has been rescaled by
  106.  *        2**(k/2) which removes the original rescaling
  107.  *        done prior to finding the mantissa square root.
  108.  *
  109.  *  NOTES
  110.  *
  111.  *    The Convergent Technologies compiler optimizes division
  112.  *    by powers of two to "arithmetic shift right" instructions.
  113.  *    This is ok for positive integers but gives different
  114.  *    results than other C compilers for negative integers.
  115.  *    For example, (-1)/2 is -1 on the CT box, 0 on every other
  116.  *    machine I tried.
  117.  *
  118.  */
  119.  
  120. /*
  121.  *  MODIFICATIONS
  122.  *
  123.  *    Function sqrt(x) is initially approximated on an
  124.  *      interval [0.5..2.0] by a quadratic polynomial with
  125.  *    the following coefficients
  126.  *
  127.  *         P0  0.3608580479718948e+00
  128.  *         P1  0.7477707028388739e+00
  129.  *         P2 -0.1105464728191369e+00
  130.  *
  131.  *    These coefficients were computed with a help of Maple
  132.  *    symbolic algebra system.  Longer approximation interval
  133.  *    is used in order to avoid multiplication by SQRT2
  134.  *    constant in case of odd exponents (this idea comes from
  135.  *    Olaf Flebbe, flebbe@tat.physik.uni-tuebingen.de).  With
  136.  *    64 bit IEEE representation three Heron's iterations are
  137.  *    good enough to satisfy essential requirements of Paranoia
  138.  *    test suite. An absolute error after three iterations is
  139.  *    below 2e-19 over the whole range (below 5e-10 after two).
  140.  *    Multiplications by powers of 2 are performed by explicit
  141.  *    calls to ldexp.
  142.  *
  143.  *    Michal Jaegermann, 21 October 1992
  144.  */
  145.  
  146. #if !defined (__M68881__) && !defined (sfp004)
  147.  
  148. #include <stdio.h>
  149. #include <math.h>
  150. #include "pml.h"
  151.  
  152. #define P0  0.3608580479718948e+00
  153. #define P1  0.7477707028388739e+00
  154. #define P2 -0.1105464728191369e+00
  155.  
  156. static char funcname[] = "sqrt";
  157.  
  158. double sqrt (x)
  159. double x;
  160. {
  161.     int k;
  162.     double m;
  163.     double u;
  164.     struct exception xcpt;
  165.     register double (*dfunc)(double, int) = ldexp;
  166.  
  167.     if (x == 0.0) {
  168.     xcpt.retval = 0.0;
  169.     } else if (x < 0.0) {
  170.     xcpt.type = DOMAIN;
  171.     xcpt.name = funcname;
  172.     xcpt.arg1 = x;
  173.     if (!matherr (&xcpt)) {
  174.         fprintf (stderr, "%s: DOMAIN error\n", funcname);
  175.         errno = EDOM;
  176.         xcpt.retval = 0.0;
  177.     }
  178.     } else {
  179.     m = frexp (x, &k);
  180.     if (k & 1) {
  181.         m = dfunc(m, 1);
  182.         /*
  183.          * multiply by 2 if our exponent was odd;
  184.          * odd k is now too big by 1 but (k >> 1) will take
  185.          * care of that
  186.          */
  187.     }
  188.     u = (P2 * m + P1) * m + P0;    /* the first guess   */
  189.  
  190.     u = dfunc((u + (m / u)), -1);    /* Heron's iteration */
  191. #ifndef __SOZOBON__
  192.     /*
  193.      * with current floating point representation used
  194.      * by Sozobon C we are already below achievable precision
  195.      */
  196.      u = dfunc((u + (m / u)), -1);   /* one more          */
  197. #endif /* __SOZOBON__ */
  198.     /*
  199.      * here we rely on the fact that -3/2 and (-3 >> 1)
  200.      * do give different results;
  201.      * the last iteration and adjust final exponent properly
  202.      */
  203.     xcpt.retval = dfunc (u + (m / u), (k >> 1) - 1);
  204.     }
  205.     return (xcpt.retval);
  206. }
  207.  
  208. double hypot(double x, double y)
  209. {
  210.     return sqrt(x*x + y*y);
  211. }
  212. #endif    /*  __M68881__, sfp004    */
  213.  
  214. #ifdef    sfp004
  215.  
  216. __asm("
  217.  
  218. comm =     -6
  219. resp =    -16
  220. zahl =      0
  221.  
  222. ");    /* end asm    */
  223.  
  224. #endif    sfp004
  225. #if defined (__M68881__) || defined (sfp004)
  226.  
  227.     __asm(".text; .even");
  228.  
  229. # ifdef    ERROR_CHECK
  230.  
  231.     __asm("
  232.  
  233. _Overflow:
  234.     .ascii \"OVERFLOW\\0\"
  235. _Domain:
  236.     .ascii \"DOMAIN\\0\"
  237. _Error_String:
  238.     .ascii \"sqrt: %s error\\n\\0\"
  239. .even
  240.  
  241. | m.ritzert 7.12.1991
  242. | ritzert@dfg.dbp.de
  243. |
  244. |    /* NAN  = {7fffffff,ffffffff}        */
  245. |    /* +Inf = {7ff00000,00000000}        */
  246. |    /* -Inf = {fff00000,00000000}        */
  247. |    /* MAX_D= {7fee42d1,30773b76}        */
  248. |    /* MIN_D= {ffee42d1,30773b76}        */
  249.  
  250. .even
  251. double_max:
  252.     .long    0x7fee42d1
  253.     .long    0x30273b76
  254. double_min:
  255.     .long    0xffee42d1
  256.     .long    0x30273b76
  257. NaN:
  258.     .long    0x7fffffff
  259.     .long    0xffffffff
  260. p_Inf:
  261.     .long    0x7ff00000
  262.     .long    0x00000000
  263. m_Inf:
  264.     .long    0xfff00000
  265.     .long    0x00000000
  266.     ");    /* end asm    */
  267. # endif    ERROR_CHECK
  268.  
  269.     __asm(".even
  270. .globl _hypot
  271. .globl _sqrt
  272. _sqrt:
  273.     ");    /* end asm    */
  274.  
  275. #endif    /* __M68881__ || sfp004    */
  276. #ifdef    __M68881__
  277.  
  278.     __asm("
  279.     fsqrtd    a7@(4), fp0    | sqrt
  280.     fmoved    fp0,a7@-    | push result
  281.     moveml    a7@+,d0-d1    | return_value
  282.     ");    /* end asm    */
  283.  
  284. #endif    __M68881__
  285. #ifdef    sfp004
  286.     __asm("
  287.     lea    0xfffa50,a0
  288.     movew    #0x5404,a0@(comm)    | specify function
  289.     cmpiw    #0x8900,a0@(resp)    | check
  290.     movel    a7@(4),a0@        | load arg_hi
  291.     movel    a7@(8),a0@        | load arg_low
  292.     movew    #0x7400,a0@(comm)    | result to d0
  293.     .long    0x0c688900, 0xfff067f8    | wait
  294.     movel    a0@,d0
  295.     movel    a0@,d1
  296.     ");    /* end asm    */
  297.  
  298. #endif    sfp004
  299. #if defined (__M68881__) || defined (sfp004)
  300. # ifdef    ERROR_CHECK
  301.     __asm("
  302. err:
  303.     lea    double_max,a0    |
  304.     swap    d0        | exponent into lower word
  305.     cmpw    a0@(16),d0    | == NaN ?
  306.     beq    error_nan    |
  307.     cmpw    a0@(24),d0    | == + Infinity ?
  308.     beq    error_plus    |
  309.     swap    d0        | result ok,
  310.     rts            | restore d0
  311. ");
  312. #ifndef    __MSHORT__
  313. __asm("
  314. error_plus:
  315.     swap    d0
  316.     moveml    d0-d1,a7@-
  317.     movel    #63,_errno    | NAN => errno = EDOM
  318.     pea    _Overflow    | for printf
  319.     bra    error_exit    |
  320. error_nan:
  321.     moveml    a0@(24),d0-d1    | result = +inf
  322.     moveml    d0-d1,a7@-
  323.     movel    #62,_errno    | NAN => errno = EDOM
  324.     pea    _Domain        | for printf
  325. ");
  326. #else    __MSHORT__
  327. __asm("
  328. error_plus:
  329.     swap    d0
  330.     moveml    d0-d1,a7@-
  331.     movew    #63,_errno    | NAN => errno = EDOM
  332.     pea    _Overflow    | for printf
  333.     bra    error_exit    |
  334. error_nan:
  335.     moveml    a0@(24),d0-d1    | result = +inf
  336.     moveml    d0-d1,a7@-
  337.     movew    #62,_errno    | NAN => errno = EDOM
  338.     pea    _Domain        | for printf
  339. ");
  340. #endif    __MSHORT__
  341. __asm("
  342. error_exit:
  343.     pea    _Error_String    |
  344.     pea    __iob+52    |
  345.     jbsr    _fprintf    |
  346.     addl    #12,a7        |
  347.     moveml    a7@+,d0-d1
  348. ");
  349. # endif    ERROR_CHECK
  350. __asm("
  351.     rts
  352.  
  353. .even
  354. _hypot:
  355.     ");
  356. #endif /* __M68881__ || sfp004    */
  357. #ifdef __M68881__
  358. __asm("
  359.     fmoved    a7@(4),fp0    |
  360.     fmulx    fp0,fp0        | x**2
  361.     fmoved    a7@(12),fp1    |
  362.     fmulx    fp1,fp1        | y**2
  363.     faddx    fp1,fp0        |
  364.     fsqrtx    fp0,fp0        | sqrt( x**2 + y**2 )
  365.     fmoved    fp0,a7@-    |
  366.     moveml    a7@+,d0-d1    | return arg
  367. ");
  368. #endif    __M68881__
  369. #ifdef    sfp004
  370. __asm("
  371.     lea    0xfffa50,a0
  372.  
  373.     movew    #0x5400,a0@(comm)    | load fp0
  374.     .long    0x0c688900, 0xfff067f8
  375.     movel    a7@(4),a0@        | load arg_hi
  376.     movel    a7@(8),a0@        | load arg_low
  377.  
  378.     movew    #0x5480,a0@(comm)    | load fp1
  379.     .long    0x0c688900, 0xfff067f8
  380.     movel    a7@(12),a0@        | load arg_hi
  381.     movel    a7@(16),a0@        | load arg_low
  382.  
  383.     movew    #0x0023,a0@(comm)
  384.     .word    0x4a68,0xfff0,0x6bfa    | test
  385.  
  386.     movew    #0x04a3,a0@(comm)
  387.     .word    0x4a68,0xfff0,0x6bfa    | test
  388.  
  389.     movew    #0x0422,a0@(comm)    | fp0 = fp0 + fp1    
  390.     .word    0x4a68,0xfff0,0x6bfa    | test
  391.  
  392.     movew    #0x0004,a0@(comm)    | sqrt(fp0)
  393.     .word    0x4a68,0xfff0,0x6bfa    | test
  394.  
  395.     movew    #0x7400,a0@(comm)    | result to d0/d1
  396.     .long    0x0c688900, 0xfff067f8
  397.     movel    a0@(zahl),d0
  398.     movel    a0@(zahl),d1
  399. ");
  400. #endif    sfp004
  401. #if ( defined (__M68881__) || defined (sfp004) ) && defined (ERROR_CHECK)
  402.  
  403. __asm("bra err");
  404.  
  405. #else    ERROR_CHECK
  406.  
  407. __asm("rts");
  408.  
  409. #endif    ERROR_CHECK
  410.